# --- Setup ---
rm(list = ls())
setwd("/Users/John/Dropbox/")

# --- Load Required Packages
library(readr)
library(interflex)
library(tidyverse)
library(ggplot2)
library(cowplot)

# --- Load Dataset ---

df <- as.data.frame(read_csv("JOP_Replication_Materials/data/processed/final_dataset.csv"))

df <- df %>% group_by(isic) %>%
  mutate(dv_2yr_lag = dplyr::lag(isic_year, 2),
         dv_3yr_lag = dplyr::lag(isic_year, 3))

# --- Interflex Model ---
intflex <- df %>%
  dplyr::select(isic, isic2, year, strategic, sqrtta, isic_year, median_share, 
                med_hhi_isic2, med_soe_isic2, dv_2yr_lag, dv_3yr_lag) %>%
  mutate(Strategic = strategic)

intflex <- as.data.frame(intflex)

# --- Poisson 2-year ---

bin.p2 <- interflex(estimator = "binning", Y = "isic_year", D = "Strategic", X = "median_share", 
                    Z = c("med_hhi_isic2", "med_soe_isic2", "dv_2yr_lag"),
                    data = intflex, vcov.type = "cluster", cl = "isic", 
                    na.rm = TRUE, main = "Marginal Effect of Strategic Status on Tech Absorption\n Poisson 2-Year Lag DV", Ylabel = "Tech Absorption Policies", 
                    Xlabel = "Median Processing Share", method = "poisson",
                    cex.main = 1, ncols=1, theme.bw = TRUE, nbins = 3, bin.labs = F)

# --- Poisson 3-year ---

bin.p3 <- interflex(estimator = "binning", Y = "isic_year", D = "Strategic", X = "median_share", 
                    Z = c("med_hhi_isic2", "med_soe_isic2", "dv_3yr_lag"),
                    data = intflex, vcov.type = "cluster", cl = "isic", 
                    na.rm = TRUE, main = "Marginal Effect of Strategic Status on Tech Absorption\n Poisson 3-Year Lag DV", Ylabel = "Tech Absorption Policies", 
                    Xlabel = "Median Processing Share", method = "poisson",
                    cex.main = 1, ncols=1, theme.bw = TRUE, nbins = 3, bin.labs = F)

# --- Linear 2-year ---

bin.l2 <- interflex(estimator = "binning", Y = "sqrtta", D = "Strategic", X = "median_share", 
                    Z = c("med_hhi_isic2", "med_soe_isic2", "dv_2yr_lag"),
                    data = intflex, vcov.type = "cluster", cl = "isic", 
                    na.rm = TRUE, main = "Marginal Effect of Strategic Status on Tech Extraction\n Linear 2-Year Lag DV", Ylabel = "Tech Absorption Policies", 
                    Xlabel = "Median Processing Share", method = "linear", 
                    cex.main = 1, ncols=1, theme.bw = TRUE, nbins = 3, bin.labs = F)


# --- Linear 3-year ---

bin.l3 <- interflex(estimator = "binning", Y = "sqrtta", D = "Strategic", X = "median_share", 
                    Z = c("med_hhi_isic2", "med_soe_isic2", "dv_3yr_lag"),
                    data = intflex, vcov.type = "cluster", cl = "isic", 
                    na.rm = TRUE, main = "Marginal Effect of Strategic Status on Tech Extraction\n Linear 3-Year Lag DV", Ylabel = "Tech Absorption Policies", 
                    Xlabel = "Median Processing Share", method = "linear", 
                    cex.main = 1, ncols=1, theme.bw = TRUE, nbins = 3, bin.labs = F)

# --- Extract figures from interflex models ---
fig_bin_p2 <- bin.p2$figure  # Poisson plot
fig_bin_l2 <- bin.l2$figure  # Linear plot

# --- Combine the plots side-by-side ---
combined_figure <- plot_grid(fig_bin_p2, fig_bin_l2, labels = NULL, nrow = 1, align = "v")

# --- Save as PDF ---
plot_path <-"JOP_Replication_Materials/appendix/output/C_figure_15.pdf"
ggsave(plot_path, combined_figure, width = 14, height = 6)

browseURL(plot_path)

# --- Extract figures from interflex models ---
fig_bin_p3 <- bin.p3$figure  # Poisson plot
fig_bin_l3 <- bin.l3$figure  # Linear plot

# --- Combine the plots side-by-side ---
combined_figure <- plot_grid(fig_bin_p3, fig_bin_l3, labels = NULL, nrow = 1, align = "v")

# --- Save as PDF ---
plot_path <-"JOP_Replication_Materials/appendix/output/C_figure_16.pdf"
ggsave(plot_path, combined_figure, width = 14, height = 6)

browseURL(plot_path)

